function [hat_beta, SE, reject, reject0, info, info0] = implement_jointtest(dy, z, shares, shifters, shifters_est, V_hh, Rh, Dtheta_h, Dtheta_dx, ind, critical , ols)

%This script implements simplified version of the joint procedure for estimation and test used in simulations
%Corresponding author: Rodrigo Adao
%Date: 09/11/2024

[Nobs, Ntheta] = size(Dtheta_dx);

%Statistic estimation
if ols == 1
   %OLS version -- normalize covariance by variance of IV
   hat_beta = ( (z'*dy)/Nobs ) /( sum(z.^2)/Nobs );  
   hat_eta = dy - hat_beta*z;
else
   %Covariance version 
   hat_beta = (z'*dy)/Nobs ;
   hat_eta = dy - ( hat_beta/( sum(z.^2)/Nobs ) )*z;
end  
hat_eta0 = dy; 

%statistics from estimated parameter
shifters_h = shifters(ind==0);
Dtheta_beta = ( Dtheta_dx'*z/Nobs )' ; 
D = - Dtheta_beta*( Dtheta_h^(-1) );

%Variance with residuals computed under the null)
R0 = hat_eta0'*shares / Nobs; 
R0_shifters_h = R0(ind==0);
   
V_bb0 = (R0.^2)*(shifters.^2); 
for nt=1:Ntheta
    V_hb0(nt,1) = (Rh(nt,:).*R0_shifters_h)*(shifters_h.*shifters_est); 
end
   
hatV0_beta = V_bb0 + D*V_hh*(D') + ( D*V_hb0 + (D*V_hb0)' );
      
%Variance with  residuals computed with estimated OLS coefficient)  
R = hat_eta'*shares / Nobs; 
R_shifters_h = R(ind==0);
  
V_bb = (R.^2)*(shifters.^2); 
for nt=1:Ntheta
    V_hb(nt,1) = (Rh(nt,:).*R_shifters_h)*(shifters_h.*shifters_est);
end
   
hatV_beta = V_bb + D*V_hh*(D') + ( D*V_hb + (D*V_hb)' );

%Standard error and tests
SE  = ( hatV_beta.^(1/2)  );
SE0 = ( hatV0_beta.^(1/2) );
if ols == 1
    SE = SE/( sum(z.^2)/Nobs );
    SE0 = SE0/( sum(z.^2)/Nobs );
end

reject = abs(hat_beta/SE) > critical;
reject0 = abs(hat_beta/SE0) > critical;
   
%Informativeness statistics
info0 = ( (V_hb0')*(V_hh^(-1))*(V_hb0) )/(V_bb0);
info  = ( (V_hb' )*(V_hh^(-1))*(V_hb ) )/(V_bb );

end